JM-ECS: A hybrid method combining the J-matrix and ECS methods for scattering 

calculations 



Y. BidasyuiQ 

Departement Wiskunde-Informatica, Universiteit Antwerpen, Antwerpen, Belgium and 
Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine 

W. Vanroose[j] J. Brocckhovej^j and F. Arick^ 

Departement Wiskunde-Informatica, Universiteit Antwerpen, Antwerpen, Belgium 



o 

> 
O 



-f— > 
I 

=5 

(N 
> 

(N 

ON 
O 

o 

> 

X 



V. Vasilevsky 

Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine 
(Dated: November 22, 2010) 
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method with exterior complex scaling and an absorbing boundary condition. The wave function is 
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in the outer region. The method is validated for a one- and two-dimensional model with partial 
wave equations, and a calculation of p-shell nuclear scattering with semi-realistic interactions. 
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I. INTRODUCTION 

Numerical simulations of breakup reactions are among 
the most complicated problems in the theoretical physics 
of nuclear reactions. The main complication arises from 
the necessity of modeling the asymptotic behavior of the 
wave function in a many-body continuum. To simplify 
the asymptotic behavior, the wave function is usually 
expressed in terms of hyperspherical harmonics [TH3]. 
But, in breakup problems, the hyperpotentials decay very 
slowly and convergence with respect to the hyper-angular 
momentum is also very slow. This results in very large 
systems of equations, even for the simplest breakup prob- 
lems. 

The Complex Scaling approach [H [5] can be used to 
avoid the problems associated with the asymptotic be- 
havior of the many-body wave function. However, this 
method is suitable only to analyze resonances in the sys- 
tem and it does not provide any breakup information, 
such as cross sections, that can be compared with exper- 
iments. 

In the numerical solution of breakup problems de- 
scribed by the Schrodinger equation it is important to 
describe each of the breakup channels in a correct way. 
Once the asymptotic form in each of the channels is un- 
derstood, it is possible to formulate a set of equations 
whose solution yields both the wave function in the in- 
teraction region, and the asymptotic amplitudes in each 
breakup channel. This is the approach taken both by the 
i?-matrix [6] and the J-matrix method [7]. 
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However, for some of the breakup channels the asymp- 
totic wave function is only known at very large distance 
from the target. As a consequence very large numeri- 
cal domains are needed to cover the interaction and the 
near-field regions. In other cases, such as the three body 
breakup in the J-matrix representation, the explicit ex- 
pression of the asymptotic wave function is not known. 
The need for an explicit asymptotic wave fuction can 
be avoided by the introduction of absorbing boundary 
conditions. These enforce an "outgoing wave" bound- 
ary condition on the numerical solution, and have been 
used successfully in atomic and molecular physics [5] and 
acoustic and electromagnetic scattering problems |9j . In- 
stead of solving in a single system both the wave function 
in the interaction region and the reaction rates, these 
methods calculate the cross section as a post-processing 
step. First the equations are solved with the absorbing 
boundary conditions and then, in a second step, the am- 
plitudes are extracted from the numerical solution. 

Absorbing boundary conditions are easy to implement 
in a grid based discretization of a partial differential equa- 
tion. They have been implemented in calculations based 
on finite difference, B-splines, finite elements etc [51 UOj. 

In nuclear physics, however, one prefers to use a L 2 - 
basis, such as the oscillator eigenstates. Such basis is 
well-adapted to the fully microscopic description of the 
compound nucleus. This paper reports on the initial ef- 
forts to introduce this approach in the context of nuclear 
few-body systems. We propose a hybrid approach with 
an I? representation of the wave function in the inter- 
action region, and with a discretized grid representation 
in the outer region. Our aim is to combine the strengths 
and benefits of both approaches. 

The paper is structured as follows. In section II, af- 
ter a short review of the ECS and J-matrix method, we 
introduce the hybrid method where the wave function is 
represented in the inner region by oscillator states and in 
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the outer region by a grid. We also discuss how the ob- 
servables are extracted from the numerical representation 
of the scattered wave. In section III we present numeri- 
cal results that validate the proposed method with model 
problems for one- and two-dimensional partial wave equa- 
tions. In section IV we show results for nuclear p-shell 
scattering. We use m = 1 and ft = 1 throughout the 
paper, except where mentioned. 



II. THE HYBRID J-MATRIX WITH ECS 
METHOD 

A. Exterior Complex Scaling as an outgoing wave 
boundary condition 

Exterior Complex Scaling (ECS) was introduced by B. 
Simon tTT] and has been initially used for the calculation 
of resonance positions and widths [5]. It has also found 
widespread application by providing an absorbing bound- 
ary condition in atomic and molecular breakup problems 
with charged particles [T2J [T3] . A review is given in [5] . 
In these problems the outgoing wave is very complicated, 
and its form depends on the effective interaction of the 
outgoing particle. It is determined by the position of the 
other charged particles involved in the breakup problem. 

The derivation of ECS is based on an analytical con- 
tinuation in the complex plane. It differs from the well- 
known Complex Scaling (CS) by starting the scaling pro- 
cedure well into the asymptotic region. 

We briefly explain why it leads to outgoing wave 
boundary conditions and apply these to the Schrodinger 
equation. 



1. Enforcing outgoing wave boundary conditions 



realized by the following mixed type boundary condition 
on the solution 

u'{L) = iku{L). (3) 

Note that this boundary condition requires the explicit 
knowledge of the wave number k. 

ECS is an alternative way to enforce the same outgo- 
ing wave condition. When we analytically continue the 
equation ([I]) to complex p e C, the general solution of 
the equation, where it is homogeneous, remains a linear 
combination of the same fundamental modes as in (|2j). 
When we impose, instead of the condition Q, a homo- 
geneous Dirichlet boundary condition in the point L'eC 
that lies inside the region where the equation is homoge- 
neous, we find that u(L') = leads to B/A = exp(2ikL') 
or 

M=exp(-2fcIm(L')). (4 ) 

So if the point L' e C is chosen such that Im(L') > 
and fclm(L') 1 then \B\ <C \A\. As a result, we have, 
effectively, enforced outgoing wave boundary conditions 
in the point L' . 

The linear combination of fundamental modes with 
coefficients A and B describes the solution everywhere 
where the equation ([I]) is homogeneous. The equation is 
homogeneous between the point L and U . So the coef- 
ficient B of the solution is still much smaller than A in 
the point L and we can conclude that we also have an 
outgoing wave boundary condition in L. 

It is important to note that, in contrast to the mixed 
boundary condition rf3l) , enforcing Dirichlet boundary 
conditions in L' does not require the knowledge of k. This 
makes it possible to describe both inelastic processes and 
outgoing boundary conditions in higher dimensions. 



Consider a one-dimensional Hclmholtz equation 

(-^ 2 -k^ju{p) = f{p) (1) 

on a domain [0, L], with constant wave number k > G 
R, homogeneous Dirichlet conditions on the left bound- 
ary, and outgoing wave boundary conditions on the right 
boundary. The right hand side, f(p), is such that it is 
zero outside the interval [0, a] with a < L. 

Let us focus this discussion on the right boundary in 
L. For a < p < L the equation is a second order homoge- 
neous differential equation with constant coefficients. In 
this region the general solution can then be written as 



i{p) = Ae lkp + Be 



-ikp 



(2) 



a linear combination of two fundamental solutions. En- 
forcing outgoing wave boundary conditions in p = L 
means that coefficient B should be zero. This can be 



2. ECS contour 

We implement this by extending the interval [0, L] with 
a contour that connects L with L' . The point V is typi- 
cally chosen such that \L\ < \L'\ and k\L'\ 1. 

Such contour can be formulated as a coordinate trans- 
formation on the interval [0, Re(L')]. 



zip) 



if(p) 



if p < L. 

if L < p < Re(L'). 



(5) 



where / is an increasing function — linear, quadratic, 
. . . — with f(Re(L')) = Im(L') and f(L) = 0. 

In particular, ECS considers a linear function, and the 
coordinate transformation is therefore written as 

{p for p< L. 

L + (p - L)e w for p > L. 
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3. Application to the Schrddinger equation 

The application of these outgoing wave boundary con- 
ditions to the Schrodinger equation is straightforward. 
Consider, as an example, the two-particle model prob- 
lem for which the radial equation is 



(T + V(p) 



1(1 



2p 



:{p) = -Vji(kp), (7) 



where T = ( — 1/2) d 2 /dp 2 denotes the kinetic energy op- 
erator, k = \/2E is the wave number, I = 0,1,2,... 
the angular momentum, V the potential that only de- 
pends on the radial coordinate p, and ji(kp) the Ricatti- 
Bessel function, which is the radial part of the incom- 
ing wave. The radial solution of this equation is ip(p) = 
jl(p)+ipsc(p), a sum of the incoming wave, ji (p) , and the 
scattered wave, ip sc (p). The latter should fit the outgoing 
wave boundary conditions. 

For problems with short range potentials the 
Schrodinger equation reduces to the Hclmholtz equation 
outside the range of the potentials with a wave number k. 
We can then apply the ECS transformation in the region 
where it reduces to a Helmholtz equation. 



4- Discretization 

ECS has been implemented in finite differences, B- 
splines and spectral elements [5]. In this section we 
discretize the differential operator that appears in the 
Schrodinger equation with finite differences using the 
Shortley-Weller formula for non- uniform grids [14]. It 
enables the discretization of the operator on the complex 
contour and uses complex valued mesh widths. Such a 
mesh mesh is illustrated on Fig. [I] 




i8/ 




Re(z) L 



Re(L') 



FIG. 1: The choice of grid points for the finite difference rep- 
resentation of the Helmholtz equation on a exterior complex 
scaled (ECS) domain. The original problem was stated on 
[0, L] with outgoing wave boundary conditions in L. The do- 
main is extended with [L, L'] where a complex grid distance 
he lB is used. In L' homogeneous Dirichlet boundary condi- 
tions are enforced. 

Consider the differential operator d 2 /dz(p) 2 on the 
contour defined by the coordinate transformation ([6]). 
We define a uniform grid 



with zq = and z n — L and mesh width h = 1/n G K, 
and a second uniform grid on the complex contour 



with z. 



(Zi)n<i<n+m On [L, L'\ 

L' and complex mesh width h 7 — (L 1 



L)/m. The union of these two grids is the grid 

(*i)o<i<n+m on [0, L] U [L,L'\ (8) 

in the entire ECS domain. To approximate the second 
derivative in each grid point Zi we use 

d 2 u(zi) 2 



dz 2 hi- 1 + hi 

x f 7— "(^-l) - (7— !~ 7T J u ( z *) + TT u ( z i+l) 

\hi-i \hi-i hi J hi 



(9) 



where hi-i and hi are the left and right mesh widths re- 
spectively, which may be complex valued. The formula 
reduces to regular second order central differences when 
hi_i = hi, i.e., in the interior real region [0, L], and in the 
interior of the complex contour [L, L'], since the scaling 
function / is taken to be linear. The only exception is 
the point z n where we lose at most an order of accuracy. 
However, with ample discretization steps, the overall ac- 
curacy is anticipated to match up to second order. Note 
that a higher order discretization at the hinge is proposed 
in |15j to maintain a second order accuracy throughout 
the domain. 

The Hamiltonian of the one-dimensional model Q is 
then a (n + m) by (n + m) matrix H/ € C n+m . 

H; = T + l l±}l diag{1/z f) + diagCVC*,)), (10) 

where T is the finite difference representation of the sec- 
ond derivative on the complex contour, and diag is a di- 
agonal matrix with the potential evaluated in each grid 
point Zi. 

For two-dimensional problems, the Hamiltonian is dis- 
cretized starting from the one-dimensional discretized 
Hamiltonian using the Kronecker product 

H 2 _d = H (l ® I + I ® H h + diag(Iq 2 (z l , zj)), (11) 

where I is a (n + m) by (n + m) unit matrix, and diag(V") 
is the two body potential evaluated at the grid points. 
The matrix H2D is now a complex valued (n + m) by 
[n + m) 2 matrix. 

The numerical solution of the PDE Q is found by 
solving the linear system 



(H l -E)x=b, 



(12) 



0<i<n 



on [0, L] 



where b is a numerical representation of the right hand 
side of the equation. 

Some examples of one- and two-dimensional wave func- 
tions on the ECS domain are presented in Fig. [2] and [3j 
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FIG. 2: The real and imaginary part of the scattered wave 
for an Helmholtz problem |T| on the ECS grid. Because of the 
homogeneous Dirichlet boundary condition the wave is only 
outgoing in L. 



potential V the solution can be written as a linear com- 
bination ip sc = Ahf (kr) + Bhf (kr) , where are the in- 
and outgoing Ricatti-Hankel functions. The coefficient A 
is then 



A = W 



(V sc (p) , hj (kpj)/w(h+(kp), hj (k P ) 



(13) 



where p € [a, L] is outside the range of the potential but 
still on the real part of the ECS domain. The Wronskian 
is calculated as W(u, v) = u'v—v'u. Since only the values 
of vj sc are known on the grid points, its first derivative is 
approximated by central differences. 

In two dimensions the problem is more complicated. 
The solution inside the numerical box has not reached its 
far-field form and a projection on the asymptotic state 
is inaccurate since the single particle and double particle 
breakup wave functions live in the same regions of space, 
especial near the edges of the domain. A procedure to 
extract the observables with the help of surface integrals 
was proposed by McCurdy, Horner and Rescigno in [16j . 
The amplitude can be written as a surface integral 




o 



FIG. 3: The real part of a 2D partial wave of a three-particle 
s-wave scattering problem on an ECS domain, (491. The 



problem has a homogeneous Dirichlet boundary condition on 
the domain [0, L'] 2 and because of the ECS this leads to a 
problem on the domain [0,L] 2 with homogeneous Dirichlet 
boundary condition on the south and west boundary, and 
outgoing wave boundary conditions on the north and east 
boundary where x = L or y = L. We show the numerical 
solution of the problem discretized with finite differences. 



5. Extraction of the observables 

Once we have the numerical solution of the equation on 
an ECS domain, we need to extract physical observables 
such as cross sections. 

In one dimension the amplitude A of the outgoing wave 
can be extracted from the numerical solution with the 
help of the Wronskian. Indeed, outside the range of the 



fh,l 2 (h,f2) =r / {T kl ,h(pi)T k2 ,i 2 (p2)'VilJ sc (p 1) p2) 



-1p sc {p 1 ,p 2 )V (T kl ,h(pi)Tk 2 ,l 2 (p2))) ■ dS, 

(14) 

where Tk 1 _i 1 (pi) is the solution of the radial equation 



2p\ 



Vi - 



r kuh (p)=0. (15) 



In a similar way r k2 ^ 2 (p) fits the same equation with T 2 , 
V2, h and fcf such that kf + k% = 2E. The contour of the 
surface integral, (14 1, should lie in the region where the 
Schrodinger equation is homogeneous. In the literature 
the functions r k 1 are often denoted as 4>u,l- I n this paper 
we have chosen a different notation to avoid confusion 
with the basis functions of the oscillator representation 
that will be considered in the next section. 

The single differential cross section (SDCS) is then 



da 
dEx 



8tt 2 



1 



U2 



\m,k 2 )Y 



(16) 



where fco is the momentum of the incoming wave |16j . 

From these amplitude formulas it is also clear why we 
use exterior complex scaling rather than complex scal- 
ing. In the latter the complete domain is rotated into 
the complex plane as p — > pe and homogeneous Dirich- 
let boundary conditions are enforced at the end of the 
grid. However, after this transformation it is very hard 
to extract scattering information. In contrast, the sole 
purpose of ECS is to provide the correct outgoing bound- 
ary conditions. It leaves the solution on the real part of 
the grid unchanged and the scattering cross sections can 
be extracted. 
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B. The J-Matrix Method 

In the original J-Matrix method for quantum scatter- 
ing [71 [T7J the wave function is represented in terms of 
some J 2 -basis set that leads to a tri-diagonal structure 
of the Hamiltonian matrix for the free-particle problem 
( Jacobi shape of the matrix) . This tri-diagonal structure 
allows the use of the complete basis set without actu- 
ally working with matrices of infinite size. A review of 
the applications of the J-matrix method can be found in 

A popular L 2 -basis set for calculations with short- 
range potentials is the set of eigenstates of the radial 
harmonic oscillator. The basis states are then generalized 
Laguerre polynomials multiplied by a weight function 



x p exp 



2b 2 



L 



1+1/2 I P 



IP 



(17) 



with a normalization 



Nit = 



2i\ 



Y{i + l + 'd/2Y 



(18) 



where i = 0, 1, 2 . . . and oscillator length b 
that is related to the oscillator frequency uj. Note that 
we have incorporated the weight p of integration coming 
from the radial coordinates in the function ( 17 1 such that 



Each basis function has a classical turning point 



R 



i,l 



by/Ai + 21 + 3. 



(19) 



The oscillator basis set is complete over L 2 and the solu- 
tions of Q can always be represented as a linear combi- 
nation of all oscillator states 



MP) =^2ci,l<f>i,l(p)- 



(20) 



i=0 



This representation reduces the Schrodinger equation to 
an infinite system of linear equations for cy. 

As already mentioned above, the kinetic energy oper- 
ator T (the free particle Hamiltonian) is a tri-diagonal 
matrix in the J-matrix method. For the oscillator basis 
the non-zero elements are 



( 2j + j + i) 



__2 



rv /(i + l)(i + i + f) 



2 



for j = i, 
for j = i — 1, 

for j = i + 



(21) 



However, the potential energy matrix in this represen- 
tation will be dense. Thus for an accurate treatment of 
the problem we need to deal with an infinitely-sized dense 
Hamiltonian matrix. As long as we deal with short-range 



potentials only, this dense potential matrix can be safely 
truncated to a finite matrix. 

This relies on the fact that highly excited oscillator 
states (with large index i) oscillate rapidly between the 
origin and the corresponding classical turning point Rij. 
So, if the range of the potential is less then its matrix 
elements will average to negligibly small value because of 
annihilating oscillatory contributions. The potential ma- 
trix can then be truncated at some i = N determined by 
the desired accuracy. Beyond this point, the Hamiltonian 
matrix is approximated by the tri-diagonal (asymptotic) 
form. We therefore refer to the dense part of the Hamilto- 
nian matrix corresponding to i < N as the "interaction 
region", and to the tridiagonal part for i > N as the 
"asymptotic region" . 

In the asymptotic region, the tri-diagonal structure of 
the matrix leads to a simple three-term recurrence rela- 
tion for the oscillator expansion coefficients of the solu- 
tion: 

■I,., if, : + {J iti - E)a + A l+1 c, l+1 =0 Vi > N (22) 

Since this is a second order recurrence relation the {cj} 
can be obtained as a linear combination of two linearly 
independent fundamental solutions. In the J-matrix 
method the regular bn and the irregular nn are cho- 
sen as fundamental solutions that can be easily obtained 
from the explicit form of the matrix J, eq. (21), 



cu = hi i + t nu Vi > N. 



(23) 



These fundamental solutions of the recurrence relation 
have a corresponding coordinate-space solution 



Mp) = Bi(p) + tM(p), 



(24) 



where Bi{p) = ThLoKi&Ap)' N'i(p) = Y^L a n i,i&,i(p) 
are usually referred to as Bessel-like and Neumann-like 
respectively, due to their asymptotic behavior. In this 
context the coefficient t in this linear combination corre- 
sponds to the scattering i-matrix. 
The full solution is then 



Mp) = xi(p) + Bi(p) + tM(p), 



(25) 



where the function Xl{p) = Si^o c ii^».K/°) ^ s nonzero 
only in the interaction region. This leads to 



Hi 



tn 



i.i 



h, l +trni 



when i < N 
when i > N 



(26) 



With this form of the solution we reduce the set of un- 
knowns to {c°j, t} and the linear system to be solved has 
iV + 1 dimensions. Solving the linear system we obtain 
simultaneously the wave function of the system, {c° ; }, 
and the scattering information, t. 



C. Introduction of the JM-ECS method 

For multi-particle scattering and reaction (breakup) 
problems it is very hard to obtain an explicit form of 
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the asymptotic solutions, also in the oscillator represen- 
tation. It is then a natural approach to avoid such ex- 
plicit forms by introducing the ECS transformation in 
the asymptotic JM region by enforcing outgoing wave 
boundary conditions. 

A formal introduction of the ECS within the oscillator 
basis representation is hard for several reasons. First, ap- 
plying the ECS coordinate transformation (|6) to oscilla- 
tor functions destroys the orthogonality as we 1 as the tri- 
diagonal structure of the kinetic energy matrix. Second, 
the convergence of the oscillator basis will be strongly 
affected. Most probably, the wave function within the 
absorbing layer will be poorly reproduced even by a huge 
truncated basis set. Finally, the coordinate transforma- 
tion makes an analytical calculation of matrix elements 
difficult. 

An alternative approach is to combine the grid and the 
oscillator representation and represent the wave function 
in the asymptotic region with finite differences. This be- 
comes possible due to a specific property of the oscilla- 
tor basis. For highly excited oscillator states the oscilla- 
tor expansion coefficients are related to the values of the 
wave function on the grid of classical turning points |19j : 



h 



b\2/R n 



*z(i? n ,,)+0 



Rn 



(27) 



This allows to couple the oscillator representation of 
the wave function with the coordinate-space representa- 
tion in the high-n region. A similar argument has been 
used for building the so-called modified J-matrix method 

(MJM) ng. 



to an large index n such that the asymptotic formula for 
the expansion coefficient, (27), can be applied. 



Again, the kinetic energy operator in this hybrid rep- 
resentation is tridiagonal since it is tridiagonal both 
in finite-differences and oscillator representations. One 
should only be careful near the matching point between 
both representations so that the asymptotic formula (27) 



is a sufficiently good approximation for representing the 
solution. 

To obtain the kinetic energy in the final point of the 
oscillator representation, we use the tridiagonal kinetic 
energy formula (21). It involves a recurrence relation 
connecting the three terms c„_2, Cn-i and c„. The lat- 
ter, the coefficient c„, is unknown. Only ip(R n ) is avail- 
able. Using the asymptotic relation ( p7| , however, we 
can calculate the required matrix element as follows: 

(Jc)„_i =J„_i in _2C„_2 + J n -i in -lCn-l 
+ Jn-UnbV^/RnMRn)- 

To calculate the kinetic energy in the first point of the 
finite difference grid, the second derivative of the wave 
function has to be known. To approximate the latter with 
a finite difference formula, one needs the wave function 
in the grid points R n -i, R n an d R n + h. We again apply 
(27) to obtain ip(Rn—i) in terms of c„_i: 



f/>"(Rn) 



c n - 1 /(b y /2/R n - 1 ) - 2tP(R n ) + ^j{R n + h) 



h 2 



The coupling between both representations around the 
matching point is sketched in Fig. |4j together with the 
terms involved to determine the correct matching. 



1. A hybrid representation of the wave function 

In the hybrid JM-ECS method we represent our one- 
dimensional wave function as a vector U' in C ra+m , where 

V = (c ,ci, . . . ,c„_i, 

^{R^),i>{R n +h),...,i>{R n + {m-l)h)). [ ' 

The first n elements represent the wave function in the 
interaction region in the oscillator representation, while 
the remaining to elements represent the wave function in 
the asymptotic region on an equidistant grid that starts 
at Rn, the n-th classical turning point, and runs up to 
R n + (to — l)h with a grid distance 

h = R n -R n -i- (29) 

We assume that the matching point, which connects the 
oscillator to finite difference representation, corresponds 



Oscillator Finite Differences 







Cn-2 c n -i 


*(Rn) V(Rn + h) 


e e 


• • 




*"(S„) 



FIG. 4: This figure illustrates how the kinetic energy matrix 
elements are calculated in the last point R n -i of the oscillator 
representation and in the first point R n of the finite difference 
representation. To calculate T applied on a solution vector 
we need to translate the oscillator representation to the grid 
and vice versa. 

The structure of the matrix representation of the ki- 
netic operator around the matching point is then 
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Jn-l.n-2 Jn-~l ,n-l Jn-l,nb\/2/ R n 



l/h 2 



( : \ 

V '■ J 



(30) 



The representation of the potential operator in the hy- 
brid JM-ECS method is more complex. Similar to the 
JM method the potential matrix in the interaction re- 
gion, covered by the oscillator representation, is dense. 
Therefore, the full hybrid potential matrix will be dense 
in the interaction region and diagonal in the asymptotic 
finite-differences region. It is clear that the computa- 
tional complexity of the problem is determined by the 
size of the interaction region because of the large number 
of potential matrix elements that needs to be calculated. 
Shifting the matching point further outward will there- 
fore strongly affect the computation time in a negative 
way. Increasing the finite-differences asymptotic part will 



have almost no effect on the latter because of the diago- 
nal potential. 

The introduction of the outgoing wave boundary con- 
dition is now straightforward by extending the finite dif- 
ference grid with an ECS contour as explained in section 
II A We use m grid points to cover both the real part 
and the ECS part of the finite difference grid. 



Similar to the finite difference representation (II) 



one can easily construct a representation for a two- 
dimensional wave function. Instead of the vector form 



(28), the solution can be represented by a matrix U' £ 
CF+m)x(n+m) that has the following structure 



coo 


■ co„-i 


don 


• d()n+m 


\ 


ClO 


■ Cl„_l 


din 


• d\ n +ni 




Cn-l a 


■ Cn — In — 1 


dn— 1 n 








■ d nn —\ 


ip(R n ,R n ) 


■ *P{Rn, Rn+rn 


) 


4+io ■ 


■ 4+ln-l 


1p(Rn+l, Rn) ■ 


■ ^{Rn+1, -Rn4 


771 ) 


dn+m ■ 


■ d n -i- mn —i 


1p(Rn+m, Rn) ■ 


■ lp{Rn+m, Rn 


fm) / 



(31) 



An example of this wave function is presented in Fig. [5] 
The spatial wave function can then, depending on the 
region in the two-dimensional domain, be written as fol- 
lows. For pi, P2 < Rn one has 

n-l 

1>(fti,P2) = E ViiMPiWM, (32) 

i,j=0 

for pi < R n and k > n one writes 

n-l 

iP(p u R l ) = J2*u<l>i(pi) (33) 

i 

and for pi < R n and I > n, 

n-l 

1>{Rk,P2)=Y,Vkj<l>M (34) 

3 

and, finally, for k,l > n 

il>{Rk,R k ) = V M . (35) 



The Hamiltonian of a 2D problem is again constructed 
as a Kronecker product of the ID Hamiltonians and a 
two-body potential. The locations of non-zero matrix 
elements are the same as in the Kronecker product of 
the two two-particle potential matrices. This leads to 
a number of dense blocks distributed over a large sparse 
matrix. In this case the computational complexity is also 
determined by the size of the finite-differences part, as 
it not only extends the diagonal, but also increases the 
number of dense blocks. 



2. Extracting Scattering information from a wave function 
in a hybrid representation 

As already mentioned at the beginning of this section, 
the extraction of observables is not always straightfor- 
ward in the JM method. This mainly comes from the 
fact that in the original formulations one has to solve the 
scattering problem and define the scattering parameters 
simultaneously. 
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R k L Rs(L r ) p. 



FIG. 5: The real part of the 2D wave function of the Three- 
particle s-wave scattering problem of ( 49 1 in the hybrid ap- 



proach with 40 oscillator states and the oscillator parameter 
b = 0.8. The values of X and Y in the oscillator region are 
chosen as the corresponding classical turning points to show 
the similar behavior of the wave function in different repre- 
sentations. Compare with Fig. j3l. Note the jump in the 
wave function near the J-matrix and finite difference bound- 
ary. This jump clearly indicates the different regions in the 
representation. 



FIG. 6: There are four region in the hybrid representation. 
In region a) p\,pi < R n and the wave function is a sum over 
product of two oscillator functions. In b) and d) one one 
coordinate, respectively p\ < R n and p2 < R n , is represented 
in the oscillator state. In region c) we use finite difference 
for both coordinates. The amplitude is calculated with the 
help of a surface integral along the path that is indicated by 
the solid line a distance Rk away from the axes. Note that 
between L and Re(L') the grid is complex valued because of 
ECS. 



In the hybrid JM-ECS approach in one dimension, once 
a solution is obtained, we still need to extract the ob- 
servables. Since the wave function is represented in the 
asymptotic region with a finite-difference representation, 
we can use a standard Wronskian technique to extract 
the scattering amplitude (see section II A 5 ) . 



Obtaining the scattering information from the solution 
in two dimensions requires the surface integral of section 
|II A5| This approach still needs to be elaborated in the 
hybrid representation. We consider for this discussion a 



contour S of the surface integral (14) that is piecewise 
parallel to one of the axes. This is illustrated in Fig. [6] 
We first integrate p-i from to a value R k £ [a, L] while 
fixing pi = Rk- The point Rk is located on the finite 
difference grid in the region where the equation is homo- 
geneous. Typically Rk is chosen a few grid points before 
the end of the real part of the grid. This procedure leads 
to a line integral parallel to the pi axis. We denote it as 

The second part of the surface integral integrates p\ 
from Rk to while keeping pi equal to Rk ■ This is a line 
integral parallel to the p\ axis and is denoted as 1\. 



The amplitude, (14), now consists of two parts 



/ il ,i 2 (fci,fc 2 ) = h + h, 



(36) 



where 



h 



T-fc 1 ,; 1 (Pi) r fe2,i 2 (- R fc) 



dj> sc {p 1 ,R k ) 
dp 2 



-ipsc(pi , RkjTk! ,h (pi) d 7 "^ 2 „ ^ ~ ) dpi 



and 



T k 1 ,h(Rk)Tk 2 ,l 2 (p2) 



-1p sc (Rk, P2)Tk 2 ,l 2 (p2 



dp2 



dip sc (Rk,p 2 ) 
dpi 
dT kuh (R k ) 
dpi 



(37) 



dp 2 . 
(38) 



This is in fact a surface integral of an in-product of two 
vector-valued functions. We can reverse the bounds on 
Ii without switching the sign of the argument, since the 
surface vector also changes direction. 

It is important to note that Ii probes the 2D wave 
function ip S c(pi, P2) where it is either represented as a 
sum over oscillator states in the first coordinate and a 
finite difference in the second or as a finite difference in 
both coordinates. The integral never probes the region 
where both coordinates are represented in the oscillator 
representation (see region (a) in Fig. [6]). 

To arrive at a practical expression for the amplitude, 
we split Ii into Jq "( - ■ ■) + Jr* (■ ■ ■) into an integration 
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from to R„, the part covered by the oscillator represen- 
tation, and an integration from R n to R k , covering the 
grid. 

In order to calculate the integrals we require the so- 



lutions Tk^iiip) and Tk 2l i 2 (p) of equation (15). We solve 
these equations also in the hybrid representation. The so- 
lutions are vectors t(ki) and tiki) G C'™ + "^ representing 
the numerical solution of ^^(p) and T k2 j 2 {p), respec- 
tively, with kf + k% = E. For pi G [0, R n ] the solution is 
a finite sum 



7fci,Zi(Pl) = ^2t(ki)i(j>i(pi) 



integral and find 

n-ln— 1 „r 71 

= -t k (k 2 ) J2 ^( fc i)^fe / <t>i{px)<t>Mdpi 

n-ln-1 .R n 

+ **(**) E E / <t>Api)<l>i{Pi)dpi 

3=0 i=0 J ° 
n— 1 

i=0 i=0 



+ C(max| / 0^-1). 



(43) 



and, similarly, the 2D wave function in region (b) of Fig. 
^ equals ^ sc (pi,Rk) = Yn=o Wik<f>i(j>i)- 

The calculation of the amplitude requires the first 
order derivatives dr kuh {R k )/dpi and dip 8C (pi,R k )/dp2. 
These will be approximated by central differences from 
t(fci) and V. We define 



t{k\)k+i - t(ki)k-i 
2ft 



and similarly for t'{k2)k- We also define 



1,1 ^i+lk — ^i-li 



2ft 



and 



/,2 _ ^ik+1 — fc-1 
2ft ' 



So the first term in the integral of I\ becomes 



Rn 



(39) 



(40) 



(41) 



/ \ ru ,di/j sc {p 1: R k ) 

TkuhKPllTkzMVtik) g— 

+ipsc(pi,Rk)T kl ,h {pi) dTk2 Qp Rk ^ ^ dpt 

(n-1 \ /n-1 

1=0 / \j=0 

''n-l \ \ 

.3=0 / Vi=0 / 



dpi, 
(42) 



We have used that J Q " <fo0j = Jjj — 

The second term in the integral I\ covers the grid and 
it is approximated using Simpson's rule 



II, , Q 

Tk 1 ,l 1 {Pl) T k 2 ,h( R k)j i — ^ S c(Pl,Rk) 



dp 2 



() 



+i)sc{pi,Rk)rk x ,h{Pi)^Q p ^ r k^l 2 {Rk)j dpt 

k k 
i= n i=n 

(44) 

where u>i is the weight of integration. The weight has 
a value Wi = ft, for all i except at i = n and i = k, 
the end points of integration, where it equals Wi = ft/2. 
Combining both parts I\ and I 2 one obtains 



i=0 



i=0 



(45) 



i=0 



i=0 



where 



1 



if i < n, 



ft/2 if i = n, 

if n < i < k, 
h/2 if £ = fc. 



(46) 



The factor 2 in (45 ) comes from the surface integral ( 14 ) 



III. NUMERICAL RESULTS 



where we have used that T k2t i 2 (R k ) = t k {k2) since the 
vector t{k) represents the solution in the hybrid repre- 
sentation. In the next step we reorder the sum and the 



To illustrate and benchmark the proposed hybrid JM- 
ECS method, we consider two model problems that are 
derived from a two-body and a three-body problem. 
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When these systems are expressed in spherical coordi- 
nates around the center of mass and expanded in spheri- 
cal harmonics one typically ends up with a large system 
of coupled radial equations, known as the partial wave 
equations. These partial waves can either be coupled ID 
or coupled 2D problems. We choose to benchmark our 
methods to model problems formulated as an uncoupled 
partial wave problem. 

The two-body problem reduces to a one-dimensional 
radial problem with the form given by Q. 

Similarly a two-dimensional radial problem is derived 
from a three-body problem and is, expressed as: 



T 1 + V 1 (p 1 ) + 



h(h 



1) 



2pi 



T 2 + V 2 ( P 2) 



h(h 



1) 



2p1 



+Vl2{p\,P2) - Y ) i'sciPUpz) 



Vk 



[V 2 (p 2 ) + V\2 (pi , P2 )} <fih ,0 (Pi )jh ( fc 0P2) 

[Vi(pi) + V 12 (p 1 ,p 2 )} (pi 2 ,o(p2)lh(koPi) 



^ sc (0, P2) = Vpi > and ^ se (pi,0) = Vp 2 
with outgoing wave boundary conditions for 



Pi 



or p 2 -> 00, 



(47) 

where p\ and p 2 are two radial coordinates represent- 
ing the distances of the first and second particle to the 
center of the coordinate system. The two body poten- 
tial is V\2(pi, p 2 ). The angular momenta l\ and I2 are 
non-negative integers. 

The total wave function is the sum of an incoming 
wave and a scattered wave. If we model an impact ion- 
ization problem, for example, the incoming wave will be 
a product of the target in the ground state, ip^.o and the 
incoming wave, ji 1 (kop 2 ), of the second particle. Taking 
into account symmetrization this is 

—j= {viifl{pi)ii 2 {k p2) + fi 2 ,o(P2)jh(koPi)) 



where tpi u o(pi) is an eigenstate of the operator (T\ + 
Vl + h(h + l)/2pi) with energy E . The incoming wave 
3h(koPi) has momentum ko such that k^/2 + Eo = k 2 /2. 

Note that the method we propose is not only applicable 
to impact problems but also to other breakup situations. 



Then the right hand side in ( 47 ) can be replaced by any 



other driving term. In photo-ionization, for example, the 
right hand side is the dipole operator p applied to the 
ground state of the two particle system ip(px,p2). 



A. One dimensional phase shift 

To validate the proposed hybrid JM-ECS technique we 
first solve the one-dimensional scattering problem with 
an attractive Gaussian potential: 

/ 1 d 2 1(1+1) x 2,2 

h — - Vne ' ° 

I 2 dx 2 2x 2 2 



#c)=0, (48) 



where the parameters of the potential were chosen as 
Vq = 0.2 and ro = 2. We determine the elastic scat- 
tering phase shift as a function of the momentum of the 
incoming particle. To estimate the accuracy of the hy- 
brid results we compare the phase shift to the one ob- 
tained with the Variable Phase Approach (VPA) [20j . 
which yields the exact result with very high accuracy for 
this simple problem. In Fig. [7] we display the results for 
oscillator length b = 0.3 and increasing matching point 
iV; this corresponds to an increasing size of the interac- 
tion region, and thus an increasing size of the truncated 
oscillator basis. It is clear from this figure that the results 
strongly depend on the choice of the matching point for 
a required accuracy. In Fig. [8] we display the results for 
different values of the angular momentum /, all seen to 
be of qualitatively comparable accuracy. 




FIG. 7: s-wave phase shift for a Gaussian potential (Vo = 0.2, 
ro = 2) as a function of wave number k for the problem 
of pj| . Shown are the hybrid JM-ECS results for differ- 
ent matching points N, and the VPA (exact) result. The 
oscillator length is b — 0.3. 

To quantify the accuracy we display an error plot in 
Fig. [9] where one notices an increasing accuracy with in- 
creasing size of the oscillator basis in the interaction re- 
gion. The error is seen to decrease in terms of N in both 
the small- and large-energy region. For large energies the 
error is mainly caused by the finite-differences approxi- 
mation to the second derivative. At small energies the 
error mainly comes from an inaccurate approximate re- 
lation ( p7| which improves as the size of the oscillator 
basis increases. 

In table |T] we display the error value for different sizes 
N of the oscillator basis and for different partial waves. 
For each partial wave we consider the fc-value where the 
error is maximal 

We observe that the error of the hybrid method has 
rather intricate and oscillatory behavior in terms of the 
size, but in general decreases with increasing oscillator 
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FIG. 8: Phase shifts for a Gaussian potential (Vb = 0.2, 
ro = 2) of equation ( 48 1 and different angular momenta 
obtained with JM-ECS (squares), compared to VPA results 
(solid lines). For all results b — 0.3 and N = 50. 




FIG. 9: The absolute error in phase shift for different sizes of 
the oscillator basis (b=0.3, 1=0) for model problem of fig. [7] 
and [7] 



basis. Full convergence is however not yet obtained. In 
any case it is seen that the approach provides relatively 
stable and accurate results for different values of the an- 
gular momenta, energy ranges, and sizes of the oscillator 
basis. 



N 


|A| (xl0- :i ) 




I = 


I = 1 


I = 2 




(k = 0.17) 


{k = 0.62) 


(k = 1.0) 


10 


6.9 


4.6 


3.6 


20 


7.5 


1.7 


1.2 


50 


3.7 


0.045 


0.039 


100 


1.4 


0.02 


0.076 


200 


0.41 


0.18 


0.0079 


300 


0.10 


0.12 


0.0026 



TABLE I: Absolute error A in the phase shift for different 
sizes of the oscillator basis (b = 0.3). 



We demonstrate this on a model breakup problem with 
a short-range potential taken from [2"T] : 



2dx 2 



1 

2dy 2 



V e~ x - V e-y 
+W e~ x ~ v - 



E) * = (49) 



with Vq — 3 and Wq = 10. These interaction parameters 
yield one bound state for the attractive "one-particle" 
potential Vq with an energy Eq — —0.411. We choose 
the energy of the incident particle to be Ei nc — 0.882 
(similarly to [21]) what leads to a total energy of the 
system E = 0.471. 

Fig. [lO] pictures the single differential cross section 
(SDCS) for the breakup after impact. The escaping par- 
ticles have an energy of 0.471 to share between them. 
Since the two particles are indistinguishable the cross sec- 
tion is symmetric around 0.471/2. In the figure we also 
compare with the results of the finite difference calcula- 
tions. We see a slight difference for equal energy sharing. 
To highlight the difference we have scaled the vertical 
axis. 

The error between the finite difference calculations and 
the hybrid method is given in table [n] where we look at 
the error for a particular energy sharing. As the number 
of oscillator states increases the results converge. 

In TablepTjwe display the difference between the results 
of the hybrid method and those of a full finite-difference 
calculation in terms of the the size of the oscillator basis 
N. 



C. More realistic s-wave benchmark with Gaussian 
potential 



B. s-wave benchmark with product two-body 
potentials 

As described above, the hybrid model introduced in 
this paper is easily extended to systems with more de- 
grees of freedom, in contrast to the original JM method. 



All potentials in (49) are exponentially decaying in 
both coordinates, and such a model is not very realistic. 
In real problems the inter-particle potential depending on 
the distance between particles should decay significantly 
slower. We consider a more realistic three-dimensional 
problem, but similar to the one described above. This 
could correspond to a scattering problem with two equal, 
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0.2 



0.3 



0.4 



FIG. 10: SCDS results of s-wave scattering for model problem 
( |49[ |, for a total energy of 0.471 (solid line — finite differences, 
dashed line — JM-ECS). The parameters of the calculation 
are iV=80 and fe=0.3. The vertical axis has been rescaled to 
emphasize the difference. 



N 


|A| 


10 


0.0937 


20 


0.0836 


40 


0.0678 


60 


0.0592 


80 


0.0184 



TABLE II: Absolute error A of the hybrid method w.r.t. a 
full finite-difference approach for the SDCS results of problem 
(49 1. The total energy is 0.471, and there is equal energy 
sharing between the particles. 



lation. In these calculations we had to increase the size 
of the finite-difference grid to cover the interaction re- 
gion, which leads to an increase increasing of the oscil- 
lator parameter (see 29 1. We considered b = 1.3 in this 
calculation. 

Overall the results are comparable. The small oscilla- 
tions in the finite difference results come from small nu- 
merical reflections at the point of complex scaling. These 
can be eliminated by using a smooth complex scaling or a 
PML [9 as an absorbing boundary. In a similar way, the 
hybrid method result also shows comparable oscillations. 

Table [ED shows the errors in the SDCS as a function 
of the oscillator basis size. Comparing the error values 
with the value of SDCS we can see that the relative error 
is considerably bigger than in previous problem. Apart 
of the increased complexity of the problem, the reason 



for this is that the asymptotic relation ( 27 ) becomes less 



accurate with an increase of the oscillator parameter b. 





0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 



light, particles and a third heavy one: 



1 . 

2 



V e Hri1 - Voe~ |r2 ' 
+WV- |ri ~ r212 - E] * = (50) 



with i"i and r2 the coordinates of the two (light) par- 
ticles with respect to the third one. Here we consider 
a Gaussian form of the potential for the interaction be- 
tween the light particles to simplify the subsequent calcu- 
lations, and to obtain a faster-decaying s-wave potential. 
The s-wave projection of this equation yields 



1 d 2 

2^2 



-Wr 



1&_ 
2(V 



Vbe 



- v e~y 

e -(x+y) 2 



2xy 



E * = (51) 



In Fig. [TT]we show the SCDS for this problem, and com- 
pare with the results from a full finite difference calcu- 



FIG. 11: SCDS results of s-wave scattering for model problem 
( |51[ ) with a Gaussian two body potential (solid line — finite 
differences, dashed line — JM-ECS). The parameters of the 
calculation are N—70 and 6=1.3. 



N 


|A| 


10 


0.3730 


20 


0.0221 


30 


0.0776 


50 


0.0685 


70 


0.0153 



TABLE III: Absolute error A of the hybrid method w.r.t. a 
full finite-difference approach for the SDCS results of problem 
(511. The total energy is 0.471, and there is equal energy 
sharing between the particles. 
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IV. APPLICATION TO NUCLEAR p-SHELL 
SCATTERING 

In this section we go beyond model problems and apply 
our method in the more realistic context of a-scattering 
in p-shell nuclei [25], in particular 8 Be [S3"] . 

This is a commonly used benchmark problem in nu- 
clear cluster physics. We will perform an a — a scat- 
tering calculation including fully microscopic states for 
the interaction region, and using a semi-realistic nucleon- 
nucleon interaction. We compare the results with those 
obtained in the JM approach. 

In the asymptotic region the system decays into two 
point particles, each corresponding to an a-particle. At 
closer distances, however, the internal structure of the 
clusters becomes of importance because of the micro- 
scopic interactions and the Pauli exchanges between nu- 
cleons. 

The cluster basis states have the form 

*„ L =A{$ 1 (a)$ 2 (a)4> nL (r)}, (52) 

where A is the antisymmetrization operator, and $12 (a) 
is a translation invariant shell-model state built up of 
s-orbitals for the a-particles. The state (p n L (r) is a three- 
dimensional harmonic oscillator state for the relative mo- 
tion of the a-clusters 



4>nL (r) = <f>n+n ,LM (r) 



(53) 



M n+no , L A~ p2/2 ^f (p 2 ) Y LM (?) 



_ Jr[ _ / 2I>+1) 

9 ~ b ' JVnL -^r(n + L + 3/2) 

We use no to denote the minimal value of the shell num- 
ber allowed by the Pauli exclusion principle 



no 



(N min -L)/2 for L<N r , 







for L > N„ 



(54) 



where iV m i n is the number of oscillator quanta in that 
shell. 



N n 



A — A, normal parity states: it = (— 1) J 
A — 3, abnormal parity states: it = (— 1) 



A+l 

(55) 

In this form n = 0, 1, ... numerates all the Pauli allowed 
states for the cluster relative motion. For all further de- 
tails, we refer to [22] and [55] , 

A common nucleon-nucleon interaction used for this 
sytem is a Volkov Nl (VI) potential, which can be writ- 
ten in the form 



the same parameters as in [23] and [22] . but omit the 
Coulomb interaction between protons to simplify the ef- 
fective asymptotic a — a interaction. 

The calculations of matrix elements in the internal re- 
gion is identical as in the JM approach. The difference 
lies in the asymptotic region, where the explicit asymp- 
totic potential, expanded in the oscillator basis, is now 
replaced by the finite differences approach, as discussed 
in section |II Cj after proper matching of both regions; 
the asymptotic a — a interaction reduces to the point- 
like effective free particle kinetic energy between the two 
clusters. 

In Fig. 



12 we show the J 71 = + , 2 + and 4 + phase shifts 



for the Be system calculated with the hybrid JM-ECS 
and JM approaches. For the latter we take results as 
they have been obtained in [23 , where a modified JM 
approach was considered, providing faster convergence 
[19] ; we refer to such calculation as M JM for consistency. 
The overall results are very close to each other. The 
only noticeable discrepancies occur at very low energies. 
The difference between MJM and JM-ECS phase shifts 
in this region is presented in Fig. [13] The oscillations in 
this difference correlate with the value of the derivative of 
the wave function in the matching point between the two 
representations. The source of the error is a "reflection" 
from this point, and an improvement of the method to 
reduce this discrepancy is currently under research. 




FIG. 12: Scattering phase shifts for the a— a system for differ- 
ent values of total angular momentum calculated with MJM 
(solid lines) and hybrid JM-ECS method (closed squares). All 
calculations were made with 80 oscillator states 



DISCUSSION AND CONCLUSION 



Vij = £V fc (l + mPQ exp{-(r lJ -/a k ) 2 }, (56) 
fe=i 



where m is the Majorana exchange parameter. The 
strength parameters V k determine the repulsive short- 
range core and the long-range attraction. We take 



This article reports on initial efforts to introduce ab- 
sorbing boundary conditions in quantum scattering prob- 
lems where the wave function is represented by a finite 
sum of oscillator states. We have introduced the exte- 
rior complex scaling (ECS) boundary conditions by con- 
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2 4 6 8 10 

E, MeV 



FIG. 13: Difference between MJM and JM-ECS scattering 
phases for the 2 + state of the a — a system. 

structing a representation of the wave function that com- 
bines the oscillator states with a grid based representa- 
tion. The oscillator representation covers the inner region 
of the problem, where the main interaction occurs, while 
the grid covers the near field. The two regions are nu- 
merically matched at an interface using an asymptotic 
expression for the expansion of oscillator states. 

The far field amplitude, which gives the cross section, 
is extracted from the numerical wave function using a 
surface integral. Extra care is required in the calculation 
of this integral in this representation. 

We have numerically solved several benchmark prob- 
lems and compared our results with literature and plain 
ECS calculations. 

Our results agree with existing methods and confirm 
that the proposed method works. However, the accu- 



racy still needs to be improved. This can be done by 
including mixed potential matrix elements from differ- 
ent representations, by using a different finite difference 
grid, or by using a higher order matching condition at the 
interface between the oscillator representation and the fi- 
nite difference representation. Also a proper treatment 
of the Coulomb interaction and other effective interac- 
tions with long-range behavior must be considered in the 
future development of the method. 

The building blocks presented in this papers are the 
starting point for a fully coupled calculations similar to 
j!3| . In such a calculation the wave function is a sum over 
li, mi and Zg, m2 of 2D functions ^i imi ,i2m 2 combined 
with Yi lTHl (fii)Y/ 2m2 (f^)- The linear system is then a 
coupled system where each diagonal block is a system 
like eq. ( |47[ ). From such a representation of the wave 
function we can extract not only the single differential 
cross section, but also the triple differential cross section, 
which gives the amplitude for certain breakup directions. 
In a similar way, the method can be applied to different 
multi-channel calculations of three-cluster systems and 
different reaction processes with light nuclei. In all such 
problems asymptotic description can be greatly simpli- 
fied if we use absorbing boundary conditions. 

We conclude that it is possible to introduce a boundary 
layer that absorbs all outgoing waves within a spectral 
basis such as the oscillator representation. 
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